# TODO: Add comment
# 
# Author: broken
###############################################################################




#import library
library(raster)
library(automap)
library(gstat)
library(rgdal)

modeFile="/home/phamha/Desktop/npp.tif"

#AOT values
aotRaster=raster(modeFile)
aot=values(aotRaster)

#get x,y coordinates
corxy=coordinates(aotRaster)
x=corxy[,'x']
y=corxy[,'y']

#cell id
totalCell=length(aotRaster)
cell = c(1:totalCell)


#data frame
table=data.frame(cell,x,y,aot)
newTable=table
#trainTable=subset(table,!is.na(aot)&aot!=-9999)
#testTable=subset(table,aot>0|is.na(aot))

#trainTable=subset(table,!is.na(aot))
trainTable=subset(table,!is.na(aot)&aot!=-9999)
testTable=subset(table,is.na(aot))

coordinates(trainTable) =~ x+y
coordinates(newTable) =~ x+y
coordinates(testTable) =~ x+y

k_rs = autoKrige(aot~1,trainTable,testTable)

plot(k_rs)

testTable$aot = k_rs$krige_output["var1.pred"]

universalPMRaster=aotRaster
universalPMValue=pred
universalPMRaster[1:totalCell]=universalPMValue

#save uk result to tiff
outFile="C:/test.tif"
writeRaster(universalPMRaster,filename=outFile,format="GTiff")

variogram = autofitVariogram(aot~1,trainTable)
#plot(variogram)
writeGDAL(testTable, "C:/c.tif")

newTable$aot[newTable$aot>0|is.na(newTable$aot)] = pred

universalPMRaster=aotRaster
universalPMValue=newTable$aot
universalPMRaster[1:totalCell]=universalPMValue

newTable

